Leonardo's rule, self-similarity and wind-induced stresses in trees 
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Examining botanical trees, Leonardo da Vinci noted that the total cross-section of branches is 
conserved across branching nodes. In this Letter, it is proposed that this rule is a consequence of 
the tree skeleton having a self-similar structure and the branch diameters being adjusted to resist 
wind-induced loads. 
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Leonardo da Vinci observed in his notebooks that 11 all 
the branches of a tree at every stage of its height when 
put together are equal in thickness to the trunk" jT] , which 
means that when a mother branch of diameter d splits 
into N daughter branches of diameters d$, the following 
relation holds on average 
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where the Leonardo exponent is A = 2. Surprisingly, 
there have been very few assessments of this rule, but the 
available data indicate that the Leonardo exponent is in 
the interval 1.8 < A < 2.3 for a large number of species 
[2j-[4] . As a matter of fact, Leonardo's rule is so natural 
to the eye that it is routinely used in computer-generated 
trees [5]. Yet, alternative analyses of the branching ge- 
ometry have been proposed based on analogies with river 
networks, bronchial trees, and arterial trees [5J. 

Two different models have been proposed to explain 
Leonardo's rule: the pipe model 7J, which assumes that 
trees are a collection of identical vascular vessels con- 
necting the leaves to the roots, and the principle of elas- 
tic similarity [5] [5] , which postulates that the deflection 
of branches under their weight is proportional to their 
length. However, none of these explanations are con- 
vincing. The first because the portion of a branch cross- 
section devoted to vascular transport (i.e. the sapwood) 
may be as low as 5% in mature trees and it seems thus 
dubious that the whole tree architecture is governed by 
hydraulic constraints. The second because the postulate 
behind elastic similarity is artificial, hard to relate to any 
adaptive advantage, and, furthermore, it seems unlikely 
that trees can respond to branch deflections. 

In this Letter, an alternative explanation is offered: 
Leonardo's rule is a consequence of trees being designed 
to resist wind-induced stresses. Plants are known to re- 
spond to dynamic loading for a long time, a phenomenon 
called thigmomorphogenesis [101 111] , In that line of 
thinking, Metzger [T2] proposed in the 19th century the 
constant-stress model. This model states that the trunk 
diameter varies such that the bending stress due to wind 
remains constant along the trunk length. The constant- 
stress model has been shown to agree with observations 



|13j . however, its implications on the whole branching 
architecture has not yet been addressed (except in the 
recent study of Lopez et al. [2]). The other important 
point is that constant-stress might not be the best de- 
sign since it implies that breakage is more likely to occur 
in the trunk or in large branches where the presence of 
defects is more probable. 

To address this problem, two equivalent analytical 
models are considered: one discrete, the fractal model, 
and one continuous, the beam model, inspired from 
McMahon & Kronauer [8] with the difference that wind 
loads are considered instead of the weight. 

The fractal model (Fig. [T^l) is constructed such that 
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where Ik and dk arc the length and diameter of a branch 
at rank k (with 1 < k < K) , TV is the number of daughter 
branches at each branching node, A is Leonardo expo- 
nent, and D is the fractal (Hausdorff) dimension of the 
tree skeleton [5] . Here, the tree skeleton is supposed to be 
self-similar such that D is uniform within the structure, 
but A can depend on k. 

The fractal dimension D has never been measured di- 
rectly on real trees. However, the fractal dimension of 
the foliage surface has been measured to lie in the inter- 
val 2.2 < Dioi. < 2.8 [TS] and, except for very particular 
architectures, it can be shown that D = Z3f i. . As al- 
ready suggested by Mandelbrot [2], it can thus be safely 
assumed that 2 < D < 3 [IS]. 

The beam model (Fig. [T]d) consists of a cantilevered 
beam whose width b and thickness h taper with the curvi- 





FIG. 1. Two analytical models: (a) the fractal tree model; 
(b) the continuous tapered beam model [8]. 
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linear coordinate s (with sq < s < 1). These two models 
can be linked using the principle sketched in Fig. [T|d. 
It consists in cutting the beam to form branches of ap- 
proximately square cross-sections. The beam thickness 
is then equivalent to the branch diameters, the ratio of 
width to thickness gives the number of branches, and 
s corresponds to the branch lengths (because the dis- 
tance between a branch and the tips is proportional to 
the branch length for infinite branching) 



N k ~ b/h ~ s" 
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Consider now two different types of wind loads: either 
a continuous loading due to the wind on the branches 
with a force per unit length q(s) ~ b or an end loading 
due to the wind in the leaves with a force Q applied in 
so equivalent to q(s) — QS(s — sq). Neglecting the wind 
incident angle and using the Eulcr-Bcrnoulli beam equa- 
tion, the curvatures k(s) resulting from the continuous 
load and the end load are found to scale respectively as 
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with E the Young's modulus and / ~ bh 3 the moment 
of inertia. The expression for continuous loading is only 
valid for s > «o |17) . Since the above scalings are equiva- 
lent at leading order, the analysis will be restricted to the 
case of end loading for simplicity. The maximum bend- 
ing stress occurs at the beam surface and is a = Ekh/2 
such that 
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The probability of fracture at a given rank k can be 
modeled by a Weibull distribution [18] to take into ac- 
count size effects 



P k = 1 - exp 



(6) 



where Vk = N k li c ird 2 ./ '4 is the volume of all branches of 
rank k, <7o is the strength of the material, Vq is an ar- 
bitrary volume taken to be Vq = t^i/4 and to is the 
Weibull's modulus (typically 5 < to < 20 for wood |19|). 
It can be shown that, for a given probability of frac- 
ture, the lightest design corresponds to the constant- 
stress model. However, since V k is decreasing with k, 
this implies that the trunk and bottom branches are more 
likely to fail. As discussed in [T3], a better design is ob- 
tained when the probability of fracture Pk is constant or 
increasing with k such that the tree can regrow after a 
big storm. The equiprobability of fracture is expressed 
as 



V k ~ hbs. 
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and it corresponds to a decreasing with s as observed in 
trees [TH [2U]. When Pk is increasing algebraically with k, 
the relation ^ still holds minor logarithmic correction. 



Now, using ([5]) and the Leonardo exponent 
and the diameter are found to depend on the fractal di- 
mension D and the Weibull's modulus to 
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In the case of infinite branching (i.e. K = 00 or so = 0), 
it gives 1.93 < A < 2.21 when 2 < D < 3, for to = 10. 
In other words, Leonardo's rule is recovered by assum- 
ing that the probability of fracture due to wind-induced 
stresses is constant. Note that, in d8k,b), constant-stress 
corresponds to the limit to — > 00. Note also that the 
number N of daughter branches at each node does not 
affect the result. 

To assess the robustness of these predictions when the 
asymmetry and stochasticity of branching, as well as 
the wind incident angle are taken into account, a three- 
dimensional numerical model has been developed. Fol- 
lowing Niklas & Kerchner [5T], a tree skeleton is recur- 
sively constructed as sketched in Fig. [2^,. 

Starting with a vertical trunk of length I = /trunk = 
1, parallel to the unit vector t and normal to the unit 
vector b, two daughter branches of lengths l\ = r±l and 
h = r 2^ are constructed in the plane normal to b such 
that their tangential unit vectors ti and t% are obtained 
by rotating t with the angles 9\ and 82 around b. The 
new normal vectors bi and b2 defining the successive 
planes of branching are then obtain by rotating b with 
an angle 7 around t! and t 2 respectively. This branching 
rule is recursively applied for K ranks, with a probability 
of branching p, yielding a tree skeleton as examplified in 
Fig. [2]d. The architecture of this skeleton is parametrized 
by the six dimensionless quantities: #1, 62, 7, ri, r 2 and 
p. This skeleton is self-similar with a fractal dimension 
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Once this tree skeleton is constructed, the diameters of 
each branch can be calculated. Assuming that the wind 
velocity Uu (where u is a unit vector) is parallel to the 





FIG. 2. Numerical tree model: (a) Sketch of the angles and 
unit vectors at a branching node; (b) Example of tree skeleton 
for 9i = -15°, 02 = 30°, 7 = 120°, n = r 2 = 0.75, p = 1, 
K = 10, D = 2.41 [as given by JgJ]. 
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FIG. 3. Deterministic tree: (a) Calculated branch diameters for the skeleton shown in Fig. [2]b; (b) Normalized average distance 
from the tips as a function of the diameter; (c) Calculated Leonardo exponent for each branching node (the horizontal bars 
show the mean value of A at each rank). In (a) and (b), the beam model corresponds to the relations (|8k,b). 



ground and uniform, the wind load on the leaves located 
at the tip of the terminal branch is 

Fleaf = \ P U 2 C X S U, (10) 

where p is the air density, C\ is a drag coefficient 
which will be taken to be 1 without loss of general- 
ity and So is the surface of the leaves assumed to be 
the square of the expected terminal branch length (i.e. 

°0 I 'trunk J- 

In addition, the wind exerts also a force on each branch 
such that, if n is the unit vector normal to both the wind 
and the branch (i.e. n = t x u/||t x u||), the force exerted 
on each branch is 

Fbranch - \ P U 2 C 2 dl\\t X U || 2 (n X t), (11) 

where C 2 is another drag coefficient taken to be 1, d and 

1 are the diameter and the length of the branch, and 
||t x u|| 2 is the square of the incident angle cosine. This 
force is applied on the branch center of mass such that its 
moment at the base of the branch is simply Mbranch = 

2 / t X Fbranch . 

Now each branch transmits the forces and moments 
applied at its extremity (either originating from upper 
branches or by leaves) such that, if F top and M top are 
the sum of forces and moments at a branch end, the force 
and moment at the branch base are 

Fbasc — Fbranch ~t~ F^op; (12a) 
Mbasc = Mbranch + M top + / 1 X F top . (12b) 

The moment at the base has two components: a 

bending moment of intensity Afbcnd = || Mbasc x t|| and 
a torsional moment of intensity M tw i st = || Mbasc ■ t||. 
The corresponding maximal bending (tensile and com- 
pressive) and shear stresses are crbcnd = ^-Mbcnd/d 3 and 

"shear = ^^twist/^ 3 - 

Assuming that there is a uniform probability of frac- 
ture (1 — e" 1 ) for every rank as given by the diame- 
ter of each branch can be calculated recursively, starting 



from the tips and ending with the trunk. In doing so, 
resistance to bending and twisting has been considered 
and the wind direction has been assumed to vary with 
increments of 45°. In this calculation, the Cauchy num- 
ber, GY = pJ7 2 /(To j bcnd, appears as the dimensionless pa- 
rameter which sets the scaling of branch diameters (but 
not their relative values) such that d ~ (j™/( 3m 2 ^ trunk . 
It has been taken to be Cy = 10 -4 which corresponds 
roughly to U = 40 ms" 1 and <r ,bcnd = 20MPa [TJ. The 
other important dimensionless numbers are the relative 
surface of leaves So (which sets the total height of the 
tree assuming that leaves have always the same dimen- 
sion whatever the size of the tree) , and the ratio of bend- 
ing to shear strength cro,bend/°'o, shear taken to be equal 
to 5 as it is generally observed for wood 19J. 

The result of such a calculation is shown in Fig. [3] for 
the deterministic skeleton pictured in Fig. [2]d. To com- 
pare these results with the theoretical predictions, the 
ratio (L)/L max is used [8], where (L) is the average dis- 
tance from the branch tips considering all possible paths 
and -L max = 1/(1 —r) is the mean ground-to-tips distance 
for an infinitely branching tree. The ratio (L) / X max is 
equivalent to (s — So) for the beam model. 

As seen in Fig. |3j the beam model accurately predicts 
the branch diameters and the Leonardo exponent. It 
means that the wind incident angle and the geometric 
details of branching do not affect these scalings. Note 
that, because of the finite number of recursions, the slope 
in Fig. is not constant as already observed in [5J . 

In Fig. gj the same simulation is run except that, at 
each branching node, the angles 9i, Q 2 , 7 are randomly 
chosen with a normal distribution of means Si, 6 2 , 7 and 
standard deviation of 10°. Same is done for r*i and r 2 of 
means f\ and f 2 and standard deviation of 0.1. It results 
that the Leonardo exponent is more scattered but the 
beam model still predicts it correctly (Fig. |4]d). 

Figure |4j; shows directly the prediction of Leonardo's 
rule by comparing the total cross-sectional areas at every 
rank. Depending on the particular angles defining the 
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FIG. 4. Stochastic tree: (a) Silhouette of a tree with parameters: 9i = -20°, 6> 2 = 40°, 7 = 100°, f\ = 0.75, f 2 = 0.85, 
p — 0.98, K — 10, D — 2.85; (b) Corresponding Leonardo exponent (same legend as in Fig. wp); (c) Total cross-sectional area 
at every rank k normalized by its maximum. The bars correspond to the tree depicted in (ajfthe curve is the prediction from 
the beam model and the symbols with the error bars correspond to the average and standard deviation on 1000 realizations 
where the mean branching angles are randomly taken in the intervals —60° < 6\ < 0, < 82 < 60°, —180° < 7 < 180°. 



branching rules, this surface can vary with a standard 
deviation of about 20% but its mean remains roughly 
constant for all ranks except the last three [35]. Remark- 
ably, the variation of this mean is accurately predicted 
by the continuous beam model. 

In summary, it has been shown that the best design to 
resist wind-induced fracture in self-similar trees naturally 
yields Leonardo's rule. The only requirement is that trees 
adapt their local growth to wind loads, a well-known phe- 
nomenon called thigmomorphogenesis whose mechanism 
at the cell level is still largely unknown. Here, the rele- 
vant property of wind loads is their divergence towards 
the branch tips, either because of the leaves or because 
the surface exposed to wind diverges. Thus the static 
loads due to the weight of fruits, snow, or ice would give 
similar results. 

The validity of the present model could be assessed 
experimentally by examining how the branch diameters 
depend on the wind peak velocity. It is predicted here 
that the following relation holds: d ~ C™^ 3m_2 ^trunk, 
where Gy is proportional to U 2 . Another assessment 
would be calculate, from a real tree skeleton, the ex- 
pected branch diameters and compare them to the mea- 
sured ones. In this Letter, aeroelastic reconfiguration 
[23], branch weight [24], and non-uniform wind profiles 
[20] , have been neglected for the sake of simplicity. It has 
also been assumed that the tree skeleton is fractal, with a 
fractal dimension 2 < D < 3. Yet, the way D and other 
features of the tree skeleton depend on the wind, and on 
the environment in general, remains to be explored. 
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